function [atlas] = formTemplates(filename, imagesPerClass, CP)
% formTemplates: 
% Parameters: CP - 2*ncp is the set of control points we have initialised
% Returns:  atlas - The cell containing the atlas data
%                altas.templateAll =  10*16*16 matrix which is the set of
%                10 image templates one for each class
%                atlas.cpAll = cell of size {10,1} with each member being a
%                matrix of size Ns*2*ncp
%                atlas.alphaAll = cell of size {10,1} with each member being a
%                matrix of size Ns*2*ncp

% initialise the subimages cell to a {10,1} cell with each cell being a
% matrix of size Ns*16*16 which contain the images
subImgs = cell(10,1);
[subImgs{10}, subImgs{1}, subImgs{2}, subImgs{3}, subImgs{4}, subImgs{5}, subImgs{6}, subImgs{7}, subImgs{8}, subImgs{9}] = readImages(filename);

% initialise the atlas data structure
atlas.templateAll = zeros(size(subImgs,1), size(subImgs{1},2), size(subImgs{1},3));
atlas.cpAll = cell(10,1);
atlas.alphaAll = cell(10,1);
for i = 1:10    
    % if we only want to do it for a smaller sample of the data
    if(imagesPerClass ~= 0)
        subImgs{i} = subImgs{i}(1:imagesPerClass,:,:);
    end    
    atlas.cpAll{i} = zeros(size(subImgs{i},1), size(CP,1), size(CP,2));
    atlas.alphaAll{i} = zeros(size(subImgs{i},1), size(CP,1), size(CP,2));
end

figure; colormap gray;
for i = 1:10
    disp(['Started constructing atlas for class number', num2str(i)]);
    atlasImages = subImgs{i};
    [template optimCP optimAlpha] = atlasFormationSubstitution(atlasImages, CP);
    atlas.templateAll(i,:,:) = template;
    atlas.cpAll{i} = optimCP;
    atlas.alphaAll{i} = optimAlpha;
    subplot(2,5,i); imagesc(squeeze(atlas.templateAll(i,:,:)));
end

disp('Done constructing altas...');

end